PHOTON DISTRIBUTION FUNCTION FOR STOCKS 
WAVE FOR STIMULATTED RAMAN SCATTERING 



O. V. Man'ko and N. V. Tcherniega 

S" 

0> ! p. N. Lebedev Physical Institute, Russian Academy of Sciences, Leninsky Pr. 53, Moscow 

. 117924, Russia 

^ ■ Abstract 



> 



New time-dependent integrals of motion are found for stimulated Raman scattering. Ex- 
plicit formula for the photon-number probability distribution as a function of the laser-field 
. intensity and the medium parameters is obtained in terms of Hermite polynomials of two 

'. variables. 

o 

a^ ; 1 Introduction 

Since the discovery of stimulated Raman scattering in 1962 [|l| this phenomenon has been inten- 
! sively investigated both theoretically and experimentally [2-7]. Quantum-mechanical description 
d \ of stimulated Raman scattering can be done in the framework of different equations, namely, 



Heisenberg-Langeven equation 0, Maxwell-Bloch equation g, and Fokker-Planck equation 
Quantum-statistical properties of stimulated Raman scattering were treated using a quadratic 
Hamiltonian For studying properties of stimulated Raman scattering with taking into ac- 
count intermolecular interaction, a qubic Hamiltonian was used |TTI . 
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_cd I Nonclassical properties of stimulated Raman scattering such as squeezing and sub-Poissonian 

statistics were considered in a number of papers [12-16]. The purpose of our work is to study the 
photon distribution function for the Stocks wave in the framework of the method of linear integrals 
of motion for systems with quadratic Hamiltonians [17-19] using the results of [20-26]. General 
formulas for matrix elements of the Gaussian density operator for the multimode oscillator in the 
Fock basis were calculated explicitly in |2^. For one-mode light described by the Wigner function 



of a generic Gaussian form with five real parameters describing the quadrature means, variances, 
and covariance, the photon distribution function was obtained explicitly in terms of Hermite 
polynomials of two variables [^ . In [^ , the temperature dependence of oscillations of the photon 



distribution function for squeezed states was investigated. It was shown that oscillations of the 
photon distribution function for squeezed and correlated light were decreasing if the temperature 
increased. 

The photon distribution function for A^-mode mixed state of light described by the Wigner 
function of the generic Gaussian form was calculated explicitly in terms of Hermite polynomials 



of 2N variables in and parameters of the photon distribution function were determined 

through the dispersion matrix and mean values of quadrature components of the light. The photon 
distribution for two-mode squeezed vacuum was investigated in |^ where its dependence on four 
parameters (two squeezing parameters, the relative phase between the two oscillators, and their 
spatial orientations) was shown. In |2^, the case of generic two-mode squeezed coherent states 
was considered, and the photon distribution function for the states was expressed both through 
four-variable and two-variable Hermite polynomials dependent on two squeezing parameters, the 
relative phase between the two oscillators, their spatial orientation, and four- dimensional shift in 
the phase space of the electromagnetic-field oscillator. 

As an application of the theoretical consideration, the linear optical transformer of photon 
statistics for multimode light was suggested |2^. The transformation coefficient was obtained 
explicitly in terms of multivariable Hermite polynomials. 

In this paper, we find new integrals of motion for the process of stimulated Raman scattering. 
We describe stimulated Raman scattering with the help of a simple model of a two-dimensional 
oscillator, the photons of the Stocks mode being described by the one mode of the oscillator and 
the phonons of the medium being described by the another mode of the oscillator. The interaction 
of the photons and phonons is taken to be quadratic in creation and annihilation operators of the 
photons and phonons. We obtain the photon-phonon probability distribution function for the 
Stocks and phonon modes after the interaction of the laser field with the medium in terms of 
Hermite polynomials of four variables. The photon distribution function for the Stocks mode is 
found explicitly and expressed both in terms of Hermite polynomials of two variables with zero 
arguments and in terms of Legendre polynomials. The mean photon number and its dispersion 
in the Stocks mode are expressed as functions of the medium parameter (temperature), the laser 
frequency, and a parameter of interaction (coupling constant). 



2 Integrals of Motion 



The simplest phenomenological Hamiltonian, which can be used for the description of one-mode 
Stocks-wave excitation, can be written 0, |2^, 



(1) 



where a and ujs are the annihilation operator and the frequency of the Stocks photon, h and uj^i 
are the annihilation operator and the frequency of the phonon, ujl is the laser frequency, and k is 
the coupling constant. The laser field is considered as classical one and its frequency is determined 
by the condition 

= ^31 + ^^S- (2) 

The damping and depletion of the laser light wave are neglected. Antistocks-mode excitation is 
also neglected and excitation of only one phonon mode is taken into consideration. 

We will show that there exist time-dependent integrals of motion for the model of Stocks-wave 
excitation. Let us construct four nonhermitian operators 



a{t) 



= ae'"^^* cosh nt + ib^e-'"^'''* sinh nt ; 



a\t) = a^e-'"^^* cosh nt - ile'"^^^^ sinh Kt ; 

h{t) = he'"^^'^ cosh Kt + m'^e"''^^* sinh ; 
P it) = fotg-i'^ait (,Qsh f^t - iae^'^st gj^h _ 

If one introduces two 4-vectors A and A (t) and uses for the 4- vectors the notation 

A = (a, b, a\ b^^ ; 
A(t) = (a{t),b{t), a\t),b\tf 



(3) 



the above relations (|^) may be represented in the matrix form 

A{t) = M (t) A. 

Here the 4x4-matrix M [t) is 



M{t) 



I e^'^st cosh Kt 

— ze*"^^^* sinh /tt e"*"^^* cosh /tt 

V -ie^'^st sinh 



^g-ic<;3it sinh Kt \ 




g-io^ait cQsh y 



(4) 



One can see that two components of the vector A, which are a, 6^ operators, are transformed 
independently of the other two components. 

In view of the commutation relations between the photon creation and annihilation operators 
a, and the phonon creation and annihilation operators 6, 6^, one can check that the operators 
constructed above satisfy boson commutation relations 

\a{t),a\t)\ = \- [6(t),6t(t)] = i, 

and the operators a(t), 6(t), and their hermitian conjugates commute as 

[a(t),6(t)] = 0; [a(t),6t(t)] = 0; 

[at(t),6(t)] = 0; [at(t),6t(t)] = o. 
It can be shown, that the total time derivatives 

da{t) dd{t) 



dt 

db{t) 
dt 

of the operators (|) are equal to zero, i.e. 

dd{t) 
dt ~ 

db{t) _ 
dt ~ 



0; 
0; 



dd\t) 

dt 
db\t) 

dt 



0; 
0. 



Consequently, the time- dependent operators a(t), d^{t), b(t), and b^{t) considered in the 
Schrodinger representation are the integrals of motion (linear with respect to the photon and 
phonon creation and annihilation operators) for Stocks-mode excitation in the framework of the 
model with the Hamiltonian (|l|). The operators (|^) are equal to standard photon and phonon 
creation and annihilation operators at the initial time moment and their commutators are time- 
independent at all time moments. 

Let us introduce quadrature components of photon and phonon creation and annihilation 
operators 

a — , a + a'^ 



b-b^ ^ b + b^ 

^^ = i7r' '' = ^- 

For Stocks-mode excitation, one can write four additional integrals of motion using the properties 
of the integrals of motion ^ 

Pa{t) = Pa cosh Ht COS Ust + % COsh Kt siu UJst 

—pb sinh hit sin u^it + qb sinh Kt cos Usit ; 

Qa(t) = —pa cosh Kt COS Ust + % COsh Kt COS OJst 

—pb sinh Kt COS tuait + % sinh Kt sin cuait ; 

(5) 

Phif) = —pa sinh Kt sin ust + % sinh Kt cos Ust 

+pb cosh Kt cos uj^it + (jb cosh Kt sin u^it ; 

Qbit) = —pa sinh Kt cos ust + (ja sinh Kt sin ust 

—pb cosh cos Usit + (Jb cosh /tt cos Ui^t . 

The physical meaning of the invariants (|^) is that their eigenvalues determine the initial values 
of classical quadrature components in the phase space of mean values {pa), {Pb), {la), and {qb). 
The number of photons does not conserve in the process of stimulated Raman scattering. But 
since any function of integrals of motion is the integral of motion |T^, one can find some 
time-dependent combinations of the photon and phonon numbers, which are integrals of motion. 
Thus, the observable 

Nait) = at(t)a(t) 

= a"'"a cosh^ Kt + (b^b + l) sinh^ Kt 



+ - |a^6^ exp [—it (lus + t^si)] — a&exp [it {lus + 1^31)] | sinh 2Kt 

is the integral of motion, which has the physical meaning of the initial number of photons in the 
system state. The observable 

Nb{t) = P{t)b{t) 

= 6^6cosh^ Kt + (ji^a + l) sinh^ Kt 

+ - |a^6^exp [—it {lus + t^si)] — abexp [it (lus + t^3i)]| sinh 



is the integral of motion, which has the physical meaning of the initial number of phonons in the 
system state. The difference of the two integrals of motion 

Na{t) - Nb{t) = d^d - b^b 

is the time-independent integral of motion, which has the physical meaning of the difference of 
photon and phonon numbers in the system, which is constant of the motion for the phenomeno- 
logical Hamiltonian (|I|). Thus, for stimulated Raman scattering we have found new integrals of 
motion. 



3 Photon— Phonon Probability Distribution Function 

In this section, we obtain explicit expression for photon distribution function of the Stocks 
mode. Let us introduce the vector column constructed from quadrature components of photon 
and phonon creation and annihilation operators in the medium at the initial time moment 

Q = {Pa^PhAaAh) 

and the vector column, constructed from the integrals of motion 

l{t) = {Pa{t),P,{t)Aa{t)Ab{t)) ■ 

Then the relation between the integrals of motion (|^) and the initial quadrature components is of 
the form 

I(t)=A(t)Q, 

where the real symplectic matrix A(t) is determined by the equation 



I cosh nt cos uJst — sinh /tt sin tuist cosh /tt sin ujst sinh nt cos uJi^t \ 

— sinh Kt sin ujst cosh Kt cos tuist sinh cos ust cosh sin uji^t 

— cosh nt cos UJst — sinh nt cos u^ist cosh K,t cos ust sinh sin uji^t 

V — sinh cos ciJ^t — cosh /tt cos c^iat svahntsinust cosh cos tuist / 



(6) 



Let us introduce the dispersion matrix of quadrature components 

a(0) 

where the matrix elements are determined through the density matrix as follows 



( ^vl 


^PaPb 


'^Paqa 


^PaQb 




^PaPh 


""pI 


'^PbQa 


^Pbqb 




'^PaQa 


^laPh 




^laqb 




V ^PaQb 


^QaPb 


^QaQb 




) 



PiP] 



<iiH3 



a. 



PiQj 



TipPiPj - {Pi){Pj); 
TTpqAj - {qi){qj); 

-TrpiqjPi +Piqj) - {Pi){qj) 



(indices i and j can be equal to a and b) . The dispersion matrix at the initial time moment t can 
be expressed through the initial dispersion matrix of quadrature components of medium photons 
and phonons in the form of matrix equation 

a{t) =A-V(0)SA^S, (7) 

where the 4x4-block matrix S consists of 2x2-zero matrices and unity matrices I2, 

^ - [I2 

If the medium photons are in the ground state and phonons are in the state of thermodynamical 
equilibrium with temperature T at the initial time moment, then 

, , 1 . X 1 

^P?(0) = 2' ^9^(0) = 2' 

^P?(0) = ^coth^; cj,2(0) = icoth^; = , 

and the matrix elements of the matrix a{t) can be found in the explicit form 
(jp2 [t) = - ^cosh^ Kt + coth ^ sinh^ nt^ ; 

(t) = - ^sinh^ nt + coth ^ cosh^ ntj ; 

Cp^p^it) = - (1 + coth — J sinh2Kt sincj^t; 



ffq-z [t) = cosh^ Kt cos cugt + - siuh^ Kt coth ^ 
o"„2 (t) = - sinh^ Kt + cosh^ Kt cos^ cui^t coth — ; 



-„^,,., — - sinh^ + cosh^ cos^ cji'^t coth — 



(8) 



cTg^gJt) = ^sinh2Kt 



cosco's't (cosc(Ji3t — sincuiat) + coth ^ cosLUi^t (cosujst — sincj^t) 



sinh^ Kt (3 cosh^ /tt 
ap^q^ (t) = ^ sm 2a;5t coth — H cos Ust (cos Ust — sm ci;5t) ; 

sinh^ Kt . cosh^ Kt B . 
^ptqb ('^J = 2 ^^"^ ^^^^ ~^ 2 ^ 2 ^^^^ "^^^^ ^^^^ ~ ^^'^ ' 

^paSi (t) = ^ sinh 2Kt ^cos {ujs — cuis) t + cos ooist coth ^ (sin oost — cos cust)^ ; 

(^ptQa (t) = ^ sinh 2Kt ^cos {us — cuis) t coth ^ + cos ugt (sin cuiat — cos tuiat)^ . 

One can see that dispersions of the photon and phonon quadratures become larger in the process of 
stimulated Raman scattering. Being initially noncorrelated the quadratures became statistically- 
dependent observables, since the Hamiltonian (|l]) is quadratic one. 



After interacting with the laser field, the state of the system can be described by the Wigner 
function of the Gaussian type 



1 



deta(t) 



exp (--Qa-\t)Q 



(9) 



where matrix a{t) is determined by formulas (P). We can express the inverse matrix cr^^(t) at 
the time moment t through the initial inverse matrix a~^{0) using (|^) and the known property of 
symplectic matrices, namely, 

a-\t) = A^a"i(0)A. 



The matrix elements of the inverse dispersion matrix a ^{t) ( ]TI1| ) have the explicit form 



Pa 



4 cosh^ Kt cos^ iUst + 2 tanh ^ sinh^ nt ; 
4 tanh ^ cosh^ nt cos^ cuist + 2 sinh^ nt ; 



sinh 2K,t 



cosLUst (coscuist — sintuist) + costuist tanh — (cosco'5t — sintu^t) 



2 cosh Kt cosLUst (sin ust — cosujst) — 2 tanh — sinh Kt sin 2uJst ; 



sinh 2Kt 



sinh 2Kt 



cos ust cosLUi^t ( 1 — tanh — ) — sincjis 



( 



P 



( . P ■ ' 
I tanh — sm ujst + cos uJst 



P 



cos uji^t COS Ust tanh 1 — sin ust sin Ui^t + tanh — cos Ui^t 



P ( 

—2 sinh^ Kt sin 2^713^ + tanh — cosh^ Kt [sm. 2a;i3t — 2 cos^ c^iat 



2 f cosh^ Kt + tanh — sinh^ Kt\ 
2 I sinh^ Kt + tanh — cosh^ Kt\ 



^qlqb^^ = sinh 2/tt sin [ujs + ^^is)^ ( 1 + tanh 



P' 



(10) 



'111 



We determine the photon-phonon probability distribution function Pnm as the probability to 
obtain n photons in the Stocks mode and m phonons in the phonon mode after the interaction of 
the laser field with the medium. The function Pnm is the matrix element of the density matrix of 
the system in the Fock basis 

Pnm = (n,m I p I n,m), 

where | n, m) is eigenstate of the set of the photon and phonon number operators c^a and 6^6 

a} a I m) = n \ n,m); 
b^b \n,m) = m\n,m). 



Using the scheme of calculations developed in [20-26] for multimode coupled oscillators we arrived 
at the expression for the photon-phonon probability distribution function Pnm in terms of Hermite 
polynomials of four variables with zero arguments 



fO 0) 



[det ((T(t) + 14/2; 



nl/2 



:i2) 



nl ml 



where the matrix R is expressed through the dispersion matrix at the time moment t 

R = (I4 - 2a{t)) (I4 + 2a{t)y^ U*, 

the matrix U is 



U 



1 

V2 



( 







-2 z 
10 10 

V 10 1/ 

and the matrix I4 is four-dimensional unity matrix. 



4 Photon Probability Distribution Function 

In experiments, the photon number in the Stocks mode is usually measured. So, it is interesting 
to average the photon-phonon probability distribution function over the phonon mode and to 
obtain the probability to have n photons in the Stocks mode. 

The photon-phonon probability distribution function can be described by the Wigner function 
d^). One has to integrate the Wigner function over the variables g?,, pb in order to obtain the 
Wigner function (averaged over the phonon mode) describing the photon state 



W^ph(ga,Pa) = ^ I I W{Q)dqbdpb 



27rJdet(T(t) 



1 



exp ( -^Qo- ^(t)Q ) dqbdpb 



(13) 



where cr~~^{t) is determined by (|TOD- The Wigner function Wp\^{qa,Pa) (0) describes the photon 
mode. 

It is convinient to change the places of quadrature components. For this purpose, we introduce 
a vector 

X = PQ, 



where the matrix P is of the form 



P 



^1 ON^ 

10 

10 

Vo ly 



Then, in view of (|T0|), for the argument of exponential function in (|T3|) one has the equality 

QcT'i(t)Q = XPA^(T~^(0)APX. 



By introducing the block matrix 



the integral in (|T3[) can be rewritten in the form 

exp (-XAX) dy, y = (pb, qt)- 



One can easily see that the integral in ( ]T3D has the form of the Gaussian integral calculated in 
Appendix A. We obtain, that the Wigner function of the photon state of the Stocks mode is 
described in explicit form by the formula 



Wpl,{pa,qa) 



exp 



where 



2 y'det a{t) det d 
(Tph' = 2a-^(J + b) d-^ (b^ + c) 



(14) 



(15) 



with the matrix elements 



-1 ^ f Sll Si2 



expressed through the matrix elements of the inverse photon-phonon matrix a ^(t) in the following 
form ????? 



Sll 



-522 



"Si 



Ob It 



PbQb 



S12 — S2I 



_ -1 _ 

Paqa 



^PaPb^PbQa ~^ ^pI ^PaQb^Qaqb ~ ^Paqb^PbQa^Pbqb ~ ^ PaPb^ PbQb^ QaQb 



(16) 



Pb it 



{^lb,b) 



The averaged probability distribution function of photons in the Stocks mode can be expressed 
through the Hermite polynomials of two variables with zero arguments using the scheme developed 
in [20-27] (see general formulas in Appendix B) 



Pn 



det ( cTph + ^ 



nl 



(17) 



where the matrix R is given by the formula 

R = U+(l2 - 2aph)(l2 + 2aph)-'U* 



with 



U= ^ 



—I I 
1 1 



One can write the photon probabihty distribution function as a function of Legendre polynomials 
using the relations between Hermite and Legendre polynomials [|rU . 
The symmetric matrix R has the form 



R 



rii ri2 
ri2 



Using the relation |T9 
where the matrix 

has the matrix element 



^i?^(0,0) = (rnr22r/^i/if (0,0) 



r 
r 



r = (riir22 



'1/2 



ri2, 



we arrive at 



det (Tph + 



-1/2 



where L„ are Legendre polynomials. 

The mean number of photons in the Stocks mode is the function of temperature and the 
coupling constant (of the laser field with the Stocks mode) and it can be calculated in explicit 
form 



in) 



cosh Kt 



COS ust ] + sinh nt coth 



(19) 



The dispersion of mean photon number in the Stocks mode is 



(T„2(t) 



1 



cosh Kt 



4 



+ COS"^ Ust + COS^ LJst 



sin LUs't 



(3 



+ - sinh nt coth — ( 1 + sin 2uJst 



/3 



-sinh sin 2ci;s'i(: coth 
2 2 



1,2 1 /5 /I 2 A 

+- smh zKt coth 77 1 - + cos Ust I 
8 2 V2 / 

1 12 . . N 1 

— - cosh Kt COS ujst (cos ujst — sm ujstj — - . 



(20) 



Thus we have found evolution of the quadrature dispersion matrix due to photon-phonon 
interaction. The propagator of the system under consideration in explicit form is presented in 
Appendix C. 



5 Conclusion 



We have shown that in the framework of simple quadratic model there exist new time- 
dependent integrals of motion for the process of stimulated Raman scattering. The linear (in 
photon and phonon quadratures) integrals of motion describe the initial values of the mean 
quadratures of the system trajectory in the phase space. The quadratic (in photon and phonon 
creation and annihilation operators) integrals of motion describe the initial numbers of photons 
and phonons in the system state. 

Another result of our work is the calculated photon distribution function, which may be ex- 
pressed either in terms of Hermite polynomials of two variables or in terms of Legendre polyno- 
mials. 

The dependence of the photon number distribution function on the parameters of the laser 
field and the couphng constant shows the possibility of processing the statistics of stimulated 
Raman scattering by varying the laser field and medium parameters (e.g., temperature). So, the 
stimulated Raman scattering can be used for production of nonclassical light. 

Analogous method of investigation can be applied for studying the stimulated Brillouin scat- 
tering. 
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Appendix A 

Here we calculate the Gaussian integral used in Section 4. We consider the integral in ([l3| ) 



02-18623. 



exp [— XAX] dyi ■ ■ ■ dy. 



where 



X=(: 



X, y) 



n 



\ym/ 



and the matrix A has four blocks 




b 
d 



We introduce notation 



k=ll=l s=lp=l fc=ls=l 



The matrix is transposed matrix c. 

Using the general formula for A^-dimensional Gaussian integral (see, for example, ) 



VdetM 

J exp i-XAX) dy = exp (-x^^x) , (21) 



J ■■■ J exp (-ZMZ + KZ) dzi-- - dz^ = ^^^^ ^ exp (^-KiVf-^K 

we get 

Vdet d 
where the matrix g is 

g = a-^(J + b) d-^ (f + c) 

and 

dy = dyi ■■ ■ di/m. 
Using the result obtained one can formulate the following theorem. 

Given an arbitrary quantum state of a composed system with a Gaussian Wigner function 
W^(Q]^,Q2) depending on the set of 2N phase-space variables Q^^ describing the first subsystem 
and on the set of 2M phase-space variables Q2 describing the second subsystem. Then the Wigner 
function of the first subsystem 

W,iQ,) = {WiQ„Q2))2, 

which is the given Wigner function of the composed system averaged over variables of the second 
subsystem Q2, has the Gaussian form. 

The statement takes place for both finite and infinite numbers of the degrees of freedom 
and M of the first and second subsystems. 



Appendix B 



In Appendix B, we derive general relations of the distribution function for a system in a 
Gaussian state to the distribution function of its subsystem. The most general mixed squeezed 
state of the A^-mode system with a Gaussian density operator g is described by the Wigner 



function W^(p, q) of the generic Gaussian form (see, for example, |T9|, |23| 



W{p,q) = (det o-)"^/^exp 



(Q))a-HQ-(Q)) 



(22) 



where 2A^-dimensional vector Q = (p, q) consists of A^ components pi, . . . ,pn and A^ components 
qi, . . . ,qN, operators p and q being the quadrature components of the creation and annihilation 
a operators (we use dimensionless variables and assume ^ = 1): 



1V2 ' 



q 



a + a' 



V2 



The 2N parameters {pi) and (g^), z = 1, 2, . . . , A^, combined into the vector (Q), are the average 
values of the quadratures, 

(p) = Tr ; (q) = Tr ^q, 
A real symmetric dispersion matrix a consists of {2N'^ + A^) variances and covariances 



a,/3 = l,2,...,2iV. 



0"a/3 — - {QaQp + QpQc 



Qa) \ Qp 



They obey certain constraints, which are nothing but the generalized uncertainty relations 
From the theorem formulated in Appendix A, it follows that the Wigner function of a subsystem 
of the A^-mode system with m degrees of freedom (1 < m < N) has the same form (p2D, namely. 



iy(p,q) = (det a)"^/2exp 



--(Q-(Q))a-i(Q-(Q)) 



(23) 



where the 2Mx2M-matrix a is obtained from the matrix a by crossing out the rows and columns 
with indices larger than m. The 2m-vector Q is obtained from the vector Q by crossing out 
components of vectors p and q with indices larger than m. It means that 



Vr(p, q) = (27r)^ Jw{p, q) dq^+i ■ ■ ■ dqN dpm+i ■ ■ ■ dp^ ■ 



(24) 



The number distribution function is nothing but the probability to have rii quanta in the first 
mode, n2 quanta in the second mode, and so on. We shall designate it by the symbol Vn, where 
the vector n consists of N nonnegative integers: n = (rii, ^2, ...,11^). This probability is given by 
the formula 

Pn = Tr ^ I n)(n |, 

where g is the density operator of the system under study, and the state | n) is a common eigenstate 



of the set of number operators dl d. 



1,2 N: 



aj I n) = nj I n). 

We use designation 

n! = ni\n2\ ■ ■ ■ un^- 
and introduce the 2A^-dimensional unitary matrix 



V2 



-Hn Hn 



The symmetric 2iV-dimensional matrix R and the 2A^-dimensional vector y, given by the relations 

R = Ut(l2w-2(T)(l2iv + 2a)-'U*; (25) 
y = 2V\l2N-2a)-\Q), (26) 

determine the number distribution. 

The factor Vq is nothing but the probability to register no quanta. It equals p3| 



det 



1-1/2 



exp 



-(Q) (2a + l2^)-^ (Q) 



If 

(Q) = o, 

the probability to have no quanta depends on {2N — 1) parameters, which coincide up to numerical 
factors with the coefficients of the characteristic polynomial of the dispersion matrix. 

The photon distribution function Vn can be expressed through the "diagonal" multidimen- 
sional Hermite polynomials 



Vn = Vo- 



H 



m 

nn 



(y) 



n! 



One can derive that 



and 



a = i (l2^ - URU*) (hN + URU*) ' 



(Q) 



:i - 2a) U 



y- 



The number distribution function for the subsystem has the same form 
ment 



Vu = Vo 



nn 



(y) 



n! 



n = (m, n2, . . . , Um) ■ 



Here 



det ( 0- + - l2m 



1-1/2 



exp 



-(Q) (2 



a 



I 



2m 



' (Q) 



(27) 
(28) 
(29) 

but with replace- 
(30) 

(31) 



The distribution function for a subsystem is completely determined by the matrix a and the vector 
Q since the matrix a and the vector Q are determined by these quantities. 

In view of the physical meaning of the joint distribution functions, one has the relation 



^n= E ■■■ E ^n. 

nm+i=0 n]v=0 

Due to this, a sum rule for the multivariable Hermite polynomials appears, namely, 

Xy) 



(32) 



•y^ ni,...,nM,n-i,...,nN^ 



(33) 



Thus, given symmetric 2A^x2A^-matrix R and the 2 A^- vector y. It means, in view of (^Sj) and 
([29|), that we have the 2A^x 2 A^- matrix a and the 2 A^- vector (Q). Using the rule formulated for 
obtaining the 2mx2m-matrix a and the 2m- vector (Q) from the matrix a and the vector (Q), we 
construct the 2m x 2m- matrix R and the 2m- vector y : 



R = VHhm-2a){hm + 2a)-'V*- 
y = 2U*(l2^-2a)-i(Q). 



(34) 
(35) 



Here U is the 2m x 2m- matrix. 



In the left-hand side of the sum rule (|33|), the constructed arguments and parameters are used. 
By multiplying the both sides of the sum rule (^) by the factor 

ni\n2\ ■■ - nj. 

and doing summation over indices rii, we arrive at the relation for generating functions for quanta 
distributions of the system and the subsystem: 



G (Ai, . . . , Xn) 



G (Ai, . . . , Am) 



OO OO \"1 \"2 

[det (AS^R + I2N 

00 00 



- M 

I ni,n2,...,npf,ni,n2 



-1/2 



00 ^rii ^712 



X"^ 11 

■^N zt{R} 



exp 
A 



...,nJR-V) 



-y(AS3 + l27v)"'S,Ay 



(36) 



ni=0 



1/2 



[det (AS^.R + hr, 



m /t{R} 

^ y(AE,R + l2, 



exp 



I ^^711, 712, .■■,flm,ni, ".2, 
1 

2 



(R^'y) 



S.Ay 



(37) 



u^u* = u*u, 



Here y = (yi, 1/2, • • • , l/2Af), the 2A^x2A^-matrix E^. is the 2A^-dimensional analog of the Pauli 
matrix ax, 

and the diagonal 2A^x2A^-matrix A reads 

N 

a = EMj' 

where each matrix A,- has only two nonzero elements: 



(A, 



jjjj 



(A, 



■j)j+N,j+N 



1. 



The matrices and A^jare defined by analogous formulas but in 2m-dimensional space. 
The relation for the generating functions is 



1) 



) Am; 1; 1; 



(3^ 



We have shown that for any Gaussian state (a state with a Gaussian Wigner function) of 
A^-mode system, the number distribution function of a m-mode subsystem is described by the 
diagonal Hermite polynomials with 2m indices. The parameters of these polynomials are deter- 
mined explicitly by quadrature dispersion matrix and quadrature means of the iV-mode system. 
The discussed model of the stimulated Raman scattering is the two-mode example of the suggested 
construction for its partial case corresponding to zero values of the quadrature means. 



Appendix C 



In this Appendix, we present the propagator of the two-mode system containing photons and 
phonons in coordinate representation. We introduce four 2x2-matrices 

^ _ / cosh nt cos ujst — sinh nt sin uJist \ 

^ \ — sinh Kt sin ust cosh nt cos uJi^t J ' 

^ _ / cosh sin tu^t sinh /tt cos cjiat \ 

^ V sinh cos cust cosh /tt sin cjiat / ' 



A3 
A4 



- cosh K,t COS ust — sinh nt cos ui^t \ 

- sinh Kt cos co'5t — cosh Kt cos cuiat / 
/ cosh Kt cos uJst sinh Kt sin uJi^t \ 

\ sinh Kt sin tu^t cosh Kt cos tOi^t J ' 



The matrix (^) determining the hnear integrals of motion i (t) is expressed in terms of these 
2 X 2-matrices 

' Ai A2 \ 
.A3 A4 i • 



A(t) 



Then, in view of general formula for the propagator of quadratic systems |T^, G (xi, X2, t) reads 
in coordinate representation 

G (xi, X2, t) = [det {-2-711X3)]'^^'^ exp j--^ X2A;^^A4X2 - 2x2Aj^xi + xiAiAj^xi | , 

where xi = {xu, xu) and X2 = (x2i, X22) are the initial and final quadrature components of the 
photon and phonon modes. 

Thus the Green function of the system under consideration is constructed. 
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